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Recently, a method to image molecular electronic wave functions using high harmonic genera- 
tion (HHG) was introduced by Itatani et al. [Nature 432, 876 (2004)]. We show that, while the 
tomographic reconstruction of general orbitals with arbitrary symmetry cannot be performed with 
long laser pulses, this becomes possible when extremely short pulses are used. An alternative recon- 
struction equation based on momentum matrix elements, rather than on dipole matrix elements, is 
proposed. We present simulations of the procedure for 2D model systems based on numerical solu- 
tions of the time-dependent Schrodinger equation, and present results from further post-processing 
of the reconstructed orbitals. 

PACS numbers: 33.80.Rv, 42.65.Ky 



High-order harmonic generation (HHG) stands for 
the emission of high-frequency radiation from a system 
driven by a strong laser field. The three-step picture [l[ 
explains HHG for atoms and molecules by a sequence of 
field ionization, acceleration of the free electron in the 
laser field and rccollision with the core. HHG has re- 
ceived increasing attention in the last years, since it can 
be used both as a source of coherent radiation, e.g., for 
creating attosecond pulses 0, 0, 0, [f| , and as a diagnostic 
tool to determine molecular properties such as the inter- 
nuclear distance [H, 0, @, • Recently, HHG has been em- 
ployed to image electronic orbitals [10( , a scheme known 
as molecular orbital tomography. There, it was shown 
that the electronic orbital of N 2 (including the sign of the 
wave function) can be reconstructed. The method is in 
principle also suited for observing femtosecond electron 
dynamics in chemical reactions. The potential of HHG 
as a femtosecond probe has been proven in experiments 
on vibrating SFg molecules 
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The main idea behind molecular orbital tomography 
is that the returning electron wave packet in the three- 
step picture can be approximately regarded as a plane 
wave. The transition matrix element describing the re- 
combination is then a Fourier transform of the bound 
state. The spectra from many different orientations of 
the molecules are combined to reconstruct a 2D projec- 
tion of the bound-state orbital on the plane orthogonal to 
the pulse propagation direction. The necessary informa- 
tion about the continuum wave packet is approximately 
obtained by measuring the spectrum of a reference atom 
with the same ionization potential as the molecular sys- 
tem. 

In a multi-electron system, corrections beyond the 
sing le-active electron model should be taken into account 
□J. As a consequence, a modified Dyson orbital and not 
the highest occupied molecular orbital is reconstructed 
HI, 14, IH. Recently it has been argued that the struc- 



orientation dependence of HHG from the more complex 
molecules acetylene and allene can still be understood in 
a single-active-electron approximation Both obser- 
vations support the idea of extending the tomographic 
technique to molecules other than N 2 . 

In the current work, we present an analysis of the 
method with important conclusions for the application 
of molecular orbital tomography to more complex sys- 
tems. In contrast to the original work [icj |. we use mo- 
mentum rather than dipole matrix elements [l8| (i.e., ve- 
locity form instead of length form) , as the velocity form 
is better suited to quantify the recombination step in 
HHG [l9| . We show that the reconstruction of an ar- 
bitrary orbital (without symmetry) is impossible within 
the original scheme based on many-cycle pulses. Instead, 
one must ensure that the electron wave packets return to 
the core from only one side. We propose to achieve this 
by using extremely short laser pulses. It is important to 
note that uni-directional rccollision is a requirement ad- 
ditional to the need of molecular head-versus-tail orienta- 
tion. Furthermore we show that using a priori assump- 
tions about the type of errors present in the method, to 
some extent those errors can be iteratively removed from 
the results. 

Let us consider a molecular system in a laser field 
polarized in the ^-direction and propagating in the z- 
direction. The intensity of the harmonic radiation with 
frequency u> is proportional to 



/M = c 2 (|PM| 2 + |P(^)| 2 ). 



(i) 



ture of the continuum wave packet does not vary between 
different rare-gas atoms 16( . Also it was shown that the 



The phases of the radiation are given by arg[P(w)]. 
Atomic units are used throughout this article. Here, 
P(oj) is the Fourier transformed dipole velocity 

P(w) = J p{t)e lut dt, (2a) 

P(*) = -(^P(x,y,z,t)\p\^(x,y,z,t)), (2b) 

with p = (—id x ,—id y ) denoting the electron momen- 
tum operator in the xy-plane. This means that we re- 
strict ourselves to harmonics emitted along the z-axis. 
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In the single-active-electron model, the time-dependent 
wave function %jj is split into two parts as ^>(x, y, z, t) = 
ipo(x, y, z, t) + ip c (x, y, z, t), where i[>o(x, y, z, t) is the ini- 
tial bound-state wave function and ip c (x, y, z, t) is the 
continuum wave packet. The time-dependence of the ini- 
tial state is given by i/}q(x, y, z, t) = ipo(x, y, z) exp (Upt), 
where I p is the ionization potential of the initial state. 

To show the need for unidirectional recollision, we first 
give a simple argument based on the laser-field-free re- 
combination of an electron wave packet. Introducing the 
plane-wave approximation, a continuum wave packet can 
be written as 



ip c (x,y,z,t) 



a{k)e 



ikx r 



0) 



where the a(k) are complex amplitudes. Neglecting de- 
pletion of the initial state, the momentum expectation 
value p(t) is given by 

p(t) ~ -(tlj (x,y,z,t)\p\ip c (x,y,z,t)) + c.c. (4) 

In the frequency domain we obtain for to > 



P(o>) = ~ f° a(k) J f ^ D (x,y)p e lkx dxdy 



x 5(w - I p - —)dk, 
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(5) 



with -0g D (a;,y) = J ipo(x,y, z)dz. The y-component of P 
is identically zero. This means that the velocity form, 
despite yielding good results for tomography, is not suit- 
able to predict the harmonic polarization correctly. For 
lo > /p, we obtain 

I(lo) = 2lo 2 \P x {u)\\ (6a) 
P x (uj) = -a{k{uj)) p*{uj) + a{-k{uj)) p(u), (6b) 

p(w) = JI ^ B (x, y)e- ik ^ x dxdy, (6c) 

where we have denned the wave number k(u>) = 
y/2(u) — I p ). To obtain p(u>) from known P x and a(k), 
we need to combine the two terms on the right-hand side 
of Eq. (|6b[) into one term, such that p(u>) is factored 
out. There are two cases in which this is possible: (i) 
the orbital is (un-)gerade or (ii) a(k) = for k > (or 
k < 0), i.e., the returning wave packets approach the nu- 
clei always from the same side. In the first case, we have 
p*{to) = ±p(u)) and the experimentally determined am- 
plitudes will be a linear combination of a{k) and a(—k). 
In this paper we focus on the second possibility. Assum- 
ing that recollision occurs only from x > 0, the Fourier 
transformed dipole velocity is given by 



P x (u) = a{-k(uj))p(u>). 



(7) 



As a modification to the above equation we use the mod- 
ified wave number k(u>) = y/2uj, as in [Io| . The modi- 
fication is related to the fact that high harmonics seem 



to be better described when the ionization potential is 
left out of the relation Q. This is also supported by 
numerical tests that we have performed, although for 
SFA-based simulations it was recently reported that the 
original relation should probably be used [ljl, [2(| ■ Phys- 
ically, the argument is that when describing the return- 
ing wave packet as a plane wave, we should take into 
account that at the moment of recombination its wave 
number is modified due to the kinetic energy being in- 
creased by an amount of the order of I p . This amounts to 
the WKB-approximation for the returning wave packet 
2l| . Although it is becoming experimentally possible to 
measure both the amplitude and the phase of the har- 
monics 0, HI , 2H 24 , [2^1 , in the original work [l(| only 
the intensities were measured. The phase information 
was added from a priori considerations. 

There are two unknowns in our set of equations, 
namely a(k(u>)) and i/jq® (x , y) . One of the main ideas 
of the procedure is to solve this problem by comparison 
with an atomic reference system for which i/j^\x,y) is 
known and a^ a '(k) is very similar to a(k). (With the 
superscript '(a)' we denote reference-system quantities.) 
This is approximately the case if the reference system 
and the molecule of interest have the same ionization po- 
tential [h}. 

In the experimental implementation (loj . the 
molecules are aligned along directions within the 
xy-plane, with an angle 8 between the molecular 
axis and the electric field. Harmonic generation 
is then considered for all orientations 9 and de- 
termined by the projections of the bound state 
il>^(x, y) = t/>q D (x cos6> + ysin9, — xsin9 + ycosO). 
The general Fourier transform 



l{ki,k 2 



^ D (x,y)e- l{klX+k2y) dxdy (8) 



can be inverted to get 



^o D (x,y) 



(27T) 5 



l {k 1 ,k 2 )e l ^ x+k ^dk 1 dk 2 . (9) 



We substitute k\ = kcosO and k 2 = — fcsin^ and change 
the integration variables accordingly, and use pe(oj) = 
g(k(u)) cos 9, —k{u>) sin 6) to arrive at 



dudO 



(10) 



Using Eq. (TTT)|) . one can reconstruct ipQ® from pe(uj). 

A more rigorous analysis of the tomography scheme is 
based on the Lewenstein model [26| . Using the veloc- 
ity form [l9l | . and a laser pulse that is turned on after 
t = 0, and turned off before t = T P , the dipole velocity 
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Px,e(u > 0) reads 

p*A u ) = ~ i I P dt I dt ' 

Jo Jo 

x v* g (k s (t,t') + A{t))e- lS{t ' t,)+lut 



(11) 



2n 



e + i{t- t') 



3/2 



where we have neglected the 'c.c.'-term in 
Eq. (|3|). The saddle-point momentum is 

given by k a (t, t') = - f t , A(t")dt"/(t - t'). Here 
A(t) = - f ^ dt'E(t') and E(t) is the electric field 
of the laser pulse. The matrix elements describing the 
ionization and recombination, and the semi-classical 
action are given by 



dion,e(M) 

Vg(k) = 

S(t,t') 



^ e { XlV )E{t)xe-^^%, (12a) 



(2tt) 3 / 2 ' 



— ikx 



id x i>o,e(x,y) 



dxdy 



dt 



„ [k s (t,t') + A(t")] 



(27r) 3 / 2: 
"M2 



(12b) 



+ I p (t-t'). (12c) 



For the tomographic procedure, it is required that in 
Eq. (fTTj) the recombination matrix element be factored 
out. To this end, we perform the integration over the 
recombination time t using the saddle-point approxima- 
tion. The saddle-point condition 



(k s (t s ,t r ) + A(t s )) 2 /2 = uj-I p 



(13) 



—ik(tu) 
(2tt) 3 / 2 



»(w) f P dt'Y, 



S(t B ,t') 



x e 



'27ri/ 



2(a; - Z p ) 



t s -t' 

x <£ on , e (fc s (t s ,0+^(*'),*') 



E(t s )k(uj) 
2tt 



(14) 



e + *(«.-#) 
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In this case molecular orbital tomography is possible and 
comparison with Eq. shows that a(—k(ui)) becomes a 
sum over all classical trajectories leading to emission of a 
photon with energy ui. Within the SFA, a(k) is identical 
for length- and velocity-form reconstruction. Note that 
the structure of Eq. pT|) . with pe(uj) as a prefactor of the 
integral, would not emerge if both positive and negative 
momenta were present (long-pulse case). 

To find a pulse for which the associated continuum 
wave packets approach the core from only one side, we 



calculate semi-classically the probability that an electron 
returns to the core with momentum k, as introduced in 
[27| . The return probability is the classical equivalent of 
|a(fc)| 2 . After tunneling, the electrons follow classical tra- 
jectories. We sample the ionization time and record the 
trajectories that return within a certain velocity interval. 
Each trajectory is weighted based on the tunneling prob- 
ability and wave-packet spreading from the Lewenstein 
model [jjj ]. A 3-cycle sin 2 -pulse with a carrier-envelope 
phase of 1.257T, for which the results are shown in Fig. [IJ 
is found to be a good candidate. To ensure that effec- 
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FIG. 1: (color online). Semi-classical distribution of return 
momenta on a linear scale (bars) and on a logarithmic scale 
(crosses) for the laser intensity 5 x 10 14 W/cm 2 and wave- 
length 780 nm. Electrons with an energy of 3.17 times the 
ponderomotive potential have k ~ ±2.57 a.u. In the inset the 
time-dependent electric field of the laser pulse is shown. 



has the form of an energy-conservation law. We as- 
sume that the electron wave packets approach the core 
always from the positive side, i.e., with negative momen- 
tum along the polarization direction. With Eq. (fT3|) . 
the return momentum is given by —k(ui), where k(u>) = duction 
y/2(u) — I p ). We obtain 



tively only a half-cycle of the pulse contributes to the 
spectrum, alternative approaches are conceivable. An 
example is using the ellipticity of the pulse, similar to 
the polarization gating method in attosecond pulse pro- 



For the simulations of orbital tomography, we study 
three cases: a symmetric orbital (ground state of 2D 
with the equilibrium internuclear distance of R = 2.0 a.u. 
as in Q ) , an antisymmetric orbital (the first excited state 
of the same system) and a prototype asymmetric orbital, 
for which we use the ground state of 2D (H-He) 2+ with 
fixed internuclear distance R = 2.2 a.u. We use a laser in- 
tensity of 5 x 10 14 W/cm 2 and wavelength 780 nm. Wc 
adjusted the softcore parameters a 2 of the softcore po- 
tentials such that I p equals 30.2 cV in all cases. For the 
asymmetric state the potential takes the form 



V(r) 



(r + f) 2 +a 2 



(r- f ) 2 + a 



(15) 



Although the results for the antisymmetric orbital could 
be improved by using harmonic phases as calculated 
from the two-center interference model @, we choose to 
present here the results using numerically calculated har- 
monic phases. First the harmonic spectra and phases are 
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FIG. 2: (color online). Simulations of orbital tomography. 
From top to bottom: the ground state of 2D with 
R = 2.0 a.u., its first excited state and the ground state of 
2D (H-He) 2+ with R = 2.2 a.u. From left to right: the exact 
bound-state orbital, the reconstructed orbital and the opti- 
mized orbital after post-processing. 



calculated by solving the TDSE for (i) many different ori- 
entations of the orbital and (ii) the reference system, 2D 
He+ for the H+ states and 2D Li 2 + for (H-He) 2+ , also 
with softcore parameters adjusted such that the ioniza- 
tion potential is 30.2 eV. Then for velocity-form recon- 
struction ag(k) is calculated according to 



(a) 



(16) 



(a) 

where Pf is the ionization yield of the reference system 
and Pi (9) is the ionization yield for different orientations 
of the molecule. The Fourier transforms pe(u) are ob- 
tained from Eq. ([7]) and the reconstruction is performed 
using Eq. (jlfjp . The mirror-symmetries present in the or- 
bitals have been used to simplify the reconstructions. Af- 
ter the tomographic reconstruction, post-processing pro- 
cedures can be applied to improve the result. We propose 
an error-reduction scheme based on the hybrid-input- 
output algorithm for phase retrieval of Fourier transforms 
introduced by Fienup [30I ]. The reconstructed image is 
transformed between Fourier and real space itcrativcly, 
and in both spaces filters are applied that follow from 
known properties of the orbital. One assumption is that 
the systematical error present in the reconstructed im- 
age has limited 8- and fc-dependence. Other important 
filters include the limited area around the origin where 
the wave function is allowed to be non-zero (the support 
area) and the realness of the orbital. We will report more 
extensively about the technique elsewhere. The results 



are shown in Fig. [2j In the cases studied here, the tomo- 
graphic reconstruction already yields orbitals that resem- 
ble well the exact orbitals. The error-reduction scheme 
improves further on these results. 

In summary, we have carried out an analysis of molec- 
ular orbital tomography. Using velocity-form recon- 
struction and uni-directional wave packets obtained with 
short pulses, orbitals of arbitrary symmetry can be 
reconstructed. A further improvement is possible by 
post-processing the reconstructed image. Molecular or- 
bital tomography represents an interesting and promising 
technique that bears great potential for observing fem- 
tosecond dynamics of electronic orbitals. 

We acknowledge discussions within the NSERC SRO 
network "Controlled electron re-scattering: femtosecond, 
sub-angstrom imaging of single molecules" . This work 
was partially supported by the Deutsche Forschungsgc- 
mcinschaft. 
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